function zero = sdp_solve(in)

global par_;

% Unpack parameters
model = par_.model;
mymoment = par_.mymoment;
mfreq = par_.mfreq;
estim = par_.estim;
show = par_.show;

lbar=NaN; % initialize parameters
ksi=NaN; %#ok<*NASGU> % initialize legacy code params

params_set;

alphapos = NaN;

switch model
        case 2  % CALVO 
                std_eA = in(1);
                menucost = menucost0;  %unused;
                alpha  = alpha0;  %unused
                omega = par_.omega;  %unused
        case 3  % GOLOSOV-LUCAS 
                menucost =  in(1); % fixed menu cost 
                std_eA = in(2);
                alpha  = alpha0;  %unused
                omega = par_.omega;  %unused
        case 4  % Woodford 
                menucost =  par_.menucost; % fixed menu cost 
                std_eA = par_.std_eA;
                alpha  = exp(in(1));   %info cost
                omega = par_.omega;  %unused
        case {5,7}  % linear hazard, uniform random mc 
                menucost =  menucost0; % unused
                omega = par_.omega;
                std_eA = par_.std_eA;
                alpha  = exp(in(1));   %slope of the hazard                
        case 6  % linear hazard 
                menucost =  menucost0; % unused
                omega = par_.omega;
                std_eA = par_.std_eA;
                alpha  = exp(in(1));   %slope of the hazard                
                alphapos  = exp(in(2));   %slope of the hazard                
end
    
if lbar>1 || lbar<0 || rho<0 || rho>1 || menucost<0 || std_eA < 0
    zero = NaN(size(in)); % tells solver to check another point
    disp('Limits reached on parameters')
    return;
end

par_.menucost = menucost;
par_.std_eA  = std_eA;
par_.alpha  = alpha;
par_.alphapos  = alphapos;
par_.omega  = omega;

matrices_set;

opts = optimset('Display','off','TolFun',sqrt(eps),'TolX',sqrt(eps));  

switch model
    case 2  % CALVO
            zero = sdp_solve_mcstd(log(par_.std_eA));
    case 3 %GL
            zero = sdp_solve_mcstd([log(par_.menucost) log(par_.std_eA)]);
    case {4,5,6,7} %Woodford solves for menucost and std_eA; linear hazard only for std_eA
            switch model
                case 4
                    in = [log(par_.menucost) log(par_.std_eA)];
                case {5,6,7}
                    in = [log(par_.omega) log(par_.std_eA)];
            end
            if estim==1
                [out,fval,flag] = fsolve('sdp_solve_mcstd',in,opts);
                sdp_solve_mcstd(out);
                if flag~=1
                    disp('sdp_solve_mcstd did not converge in sdp_solve')
                    par_.omega = omega0;            %set parameters back to starting values
                    par_.std_eA = std_eA0;
                end
            else %this just solves used for USvol; does not check whether freq and freq_target equals
                sdp_solve_mcstd(in);
            end
            switch mymoment 
                case 'kurtosis'
                    zero = log(par_.KurtPchanges/kurt_target);
                case {'hazard','hazard_USvol'}
                    %interpolate hazard and density
                    hazard_model = interp1(-par_.Pdistans,par_.lambda_prime,deviation_data,'linear','extrap');
                    density_model = interp1(-par_.Pdistans,par_.Pdist_sh,deviation_data,'linear','extrap');
                    density_model = density_model/sum(density_model);  %make sure it integrates to 1
                    nump_data = length(deviation_data);
                    switch model
                        case {2,3,4,5,7}
                            %calculate the data-density-weighted squared hazard
                            %difference, excluding the first and the last point                  
                            zero(1,1) = sqrt(density_data(2:nump_data-1)'*(hazard_model(2:nump_data-1)-freq_data(2:nump_data-1)).^2);
                        case 6
                            nump_firstposgap_data = find(deviation_data>=0,1);
                            zero(1,1) = sqrt(density_data(2:nump_firstposgap_data-1)'/sum(density_data(2:nump_firstposgap_data-1))*...
                                            (hazard_model(2:nump_firstposgap_data-1)-freq_data(2:nump_firstposgap_data-1)).^2);
                            zero(1,2) = sqrt(density_data(nump_firstposgap_data:nump_data-1)'/sum(density_data(nump_firstposgap_data:nump_data-1))*...
                                                (hazard_model(nump_firstposgap_data:nump_data-1)-freq_data(nump_firstposgap_data:nump_data-1)).^2);
                    end
                case 'stdpcdist'
                      stand_pc_dist_model = interp1((par_.Pchangegrid-par_.EPchange)/par_.StdPchanges,par_.denspchanges,std_ref_price_data,'linear','extrap');
                      stand_pc_dist_model = stand_pc_dist_model/sum(stand_pc_dist_model);  %make sure it integrates to 1
                      nump_data = length(deviation_data);
                      %difference, excluding the first and the last point                  
                      %zero(1,1) = sqrt(std_ref_price_density_data(2:nump_data-1)'*(stand_pc_dist_model(2:nump_data-1)-std_ref_price_density_data(2:nump_data-1)).^2);
                      zero(1,1) = sqrt(sum((stand_pc_dist_model(2:nump_data-1)-std_ref_price_density_data(2:nump_data-1)).^2));
            end
            fprintf('\n\nParameter alpha tried        : %0.8g   \n',alpha);
            fprintf('\n\nParameter omega tried        : %0.8g   \n',omega);
            fprintf('Distance from target         : %0.8g   \n\n',zero);
end
end